miRNA Differential expression analysis for plasma samples

Plasma only samples

import pandas as pd
samples = pd.read_csv('data/sample_sheet.csv')
# must be Healthy Control Study and must pass MiSeq QC
samples = samples.loc[(samples['Study'] == 'Healthy Controls') & (samples['MISEQ.QC.PASS'] == 'PASS')]  
samples = samples.set_index('MT.Unique.ID').sort_values(by=['Participant.ID', 'Source'])
# capitalize & strip whitespace for consistency
for column in ['Gender', 'Race', 'Source']:
    samples[column] = samples[column].str.capitalize()
    samples[column] = samples[column].str.strip()
# use correct ontology terms
race_ontology = {'Asian': 'Asian',
 'Black or african american': 'African_American',
 'Mixed/asian & white': 'Multiracial',
 'Mixed/asian &black': 'Multiracial',
 'Mixed/black, white, asian': 'Multiracial',
 'Native hawiian or other pacific islander': 'Pacific_Islander',
 'Pacific islander': 'Pacific_Islander',
 'White': 'White'}
for id in samples.index:
    race = samples.at[id, 'Race']
    samples.at[id, 'Race'] = race_ontology[race] if race in race_ontology else 'Multiracial'
# get only the plasma samples
mir_counts = pd.read_csv("data/get_canonical/canon_mir_counts.csv", index_col=0)
samples.index = pd.Index(['X' + str(row) for row in samples.index])
plasma_samples = samples.loc[samples['Source'] == "Plasma"]
plasma_mir_counts = mir_counts[plasma_samples.index]
plasma_samples.to_csv('data/plasma_samples.csv')
plasma_mir_counts.to_csv('data/plasma_mir_counts.csv')

Load the sample data and miRNA counts

samples <- read.csv('data/plasma_samples.csv')
counts <- read.csv('data/plasma_mir_counts.csv')
samples <- subset(samples, Library.Generation.Set != "SetRecheck" )  # exclude SetRecheck samples
samples <- samples[samples$X != "X11" & samples$X != "X92",] # remove outliers
samples$Participant.ID <- factor(samples$Participant.ID)  # ID is categorical, not numerical
rownames(counts) = counts$X
rownames(samples) = samples$X
counts$X <- NULL  # remove extra column
samples$X <- NULL
counts <- counts[,unique(rownames(samples))]
head(samples)
head(counts)

Filter

library(edgeR)
design <- model.matrix(~0+Race+Gender+Age, samples)
dge = DGEList(counts = counts, samples = samples)
# require miRNAs to have CPM > 1 in at least 2 samples
countsPerMillion <- edgeR::cpm(dge)
countCheck <- countsPerMillion > 1
head(countCheck)
keep <- which(rowSums(countCheck) >= 2) 
dge <- dge[keep,]

Explore variance

library(SingleCellExperiment)
library(scater)

Attaching package: ‘scater’

The following object is masked from ‘package:EDASeq’:

    plotRLE

The following object is masked from ‘package:S4Vectors’:

    rename

The following object is masked from ‘package:limma’:

    plotMDS

The following object is masked from ‘package:stats’:

    filter

reads_sce <- SingleCellExperiment(assays=list(counts=dge$counts),  colData=dge$samples)
# remove unexpressed miRNAs
keep_feature <- rowSums(counts(reads_sce) > 0) > 0
reads_sce <- reads_sce[keep_feature, ]
reads_sce <- calculateQCMetrics(reads_sce)

Note that the names of some metrics have changed, see 'Renamed metrics' in ?calculateQCMetrics.
Old names are currently maintained for back-compatibility, but may be removed in future releases.

# log transform
cpm(reads_sce) <- calculateCPM(reads_sce)
reads_sce <- normalize(reads_sce)

using library sizes as size factors

logcounts(reads_sce) <- log2(calculateCPM(reads_sce) + 1)
# visualize
hist(reads_sce$total_counts, breaks=100)  # counts per sample

hist(reads_sce$total_features, breaks=100)  # counts per miRNA

plotQC(reads_sce, type = "highest-expression")

plotQC(reads_sce, type="explanatory-variables", variables=c("Index", "Participant.ID", "Collection.Date", "Library.Generation.Set", "MiSeq.QC.Run", 'total_features', "Age", "Race", "Gender"))

Examine sources of variation

plotPCA(reads_sce, exprs_values = "logcounts", colour_by = "Library.Generation.Set", size_by = "total_features")

non-plotting arguments like 'exprs_values' should go in 'run_args'

plotPCA(reads_sce, exprs_values = "logcounts", colour_by = "Race", shape_by="Gender", size_by = "total_features")

non-plotting arguments like 'exprs_values' should go in 'run_args'

for (var in c("total_features", "Age", "Library.Generation.Set", "Index", "Participant.ID", "Race", "Gender", "Collection.Date")) {
  print(
    plotQC(reads_sce, type = "find-pcs", exprs_values = "logcounts", variable = var)
    )  + ggtitle(var) 
  }

Normalize & Remove unwanted sources of variation

library(RUVSeq)
library(ggplot2)
library(mvoutlier)

Loading required package: sgeostat
sROC 0.1-2 loaded

dge <- calcNormFactors(dge)
dge <- estimateGLMCommonDisp(dge, design)
dge <- estimateGLMTagwiseDisp(dge, design)
fit <- glmFit(dge, design)
res <- residuals(fit, type="deviance")
set <- newSeqExpressionSet(dge$counts, phenoData=dge$samples)
ruvr_sets <- list()
for(k in 1:5) {
  ruvr_sets[[k]] <- RUVr(set, row.names(dge), k=k, res)
  assay(reads_sce, paste("RUVr k=", toString(k))) <- log2(t(t(assayData(ruvr_sets[[k]])$normalizedCounts) / colSums(assayData(ruvr_sets[[k]])$normalizedCounts) * 1e6) + 1)
}
for(n in assayNames(reads_sce)) {
  print(
        plotPCA(
            reads_sce,
            colour_by = "Library.Generation.Set",
            size_by = "total_features",
            exprs_values = n
        ) + ggtitle(paste("Batch",n))
  )
  print(
        plotPCA(
            reads_sce,
            colour_by = "Race",
            shape_by = "Gender",
            size_by = "total_features",
            exprs_values = n
        ) + ggtitle(paste("Demographics", n))        
  )
}

non-plotting arguments like 'exprs_values' should go in 'run_args'

Detect outliers

reads_sce <- runPCA(reads_sce, use_coldata = TRUE, detect_outliers = TRUE)

failed to find 'pct_counts_feature_control' in column metadatafailed to find 'total_features_by_counts_feature_control' in column metadatafailed to find 'log10_total_counts_endogenous' in column metadatafailed to find 'log10_total_counts_feature_control' in column metadata

outliers <- colnames(reads_sce)[reads_sce$outlier]
head(outliers)

character(0)

Examine sources of variance after removing unwanted variation

for (var in c("total_features", "Age", "Library.Generation.Set", "Index", "Participant.ID", "Race", "Gender", "Collection.Date")) {
  print(
    plotQC(reads_sce, type = "find-pcs", exprs_values = "RUVr k= 2", variable = var)
    )  + ggtitle(var) 
}

Visualize top highly-expressed miRNAs male-vs-female

Visualize top highly-expressed miRNAs by Race

# Rank the mean expression values for plasma/serum miRs. Highest expression = 1
mean.expr.afr.rank <- rank(-1*rowMeans(norm.expr.matr[, colData(reads_sce)$Race=="African_American"]))
mean.expr.asn.rank <- rank(-1*rowMeans(norm.expr.matr[, colData(reads_sce)$Race=="Asian"]))
mean.expr.wht.rank <- rank(-1*rowMeans(norm.expr.matr[, colData(reads_sce)$Race=="White"]))
top_N <- 20
top.miRs <- row.names(norm.expr.matr)[mean.expr.afr.rank <= top_N | mean.expr.asn.rank <= top_N | mean.expr.wht.rank <= top_N] # get the names of the top miRs in plasma or serum
norm.expr.top <- norm.expr.matr[top.miRs, ] # Get the expression matrix for the top mIRs
norm.expr.melt <- reshape2::melt(norm.expr.top) # Convert your normalized expression matrix to a 3 column data.frame (row.name, col.name, expression value). Melt is in the dpylr package, I believe.
colnames(norm.expr.melt) <- c("miR.ID", "MT.Unique.ID", "norm.expr") # just so it's easier for me to tell you which columns I'm using.
norm.expr.melt$Race <- ""
for (row_num in 1:nrow(norm.expr.melt)){  # Pull the Gender values from the column metadata.
  mt_unique_id <- norm.expr.melt[row_num, ]$MT.Unique.ID
  norm.expr.melt[row_num, "Race"] <- as.character(samples[mt_unique_id,"Race"])
}
norm.expr.melt <- norm.expr.melt[norm.expr.melt$Race %in% c("African_American", "Asian", "White"),]
# This would be a simple boxplot with plasma/serum side-by-side. Overlaying the boxes over points takes a little more tweaking to get the dodge/width right, but it's doable.
ggplot(norm.expr.melt, aes(x=reorder(miR.ID, norm.expr, FUN="median"), y=norm.expr, fill=Race)) + geom_boxplot(pos="dodge", outlier.size=0.5) + 
  ggtitle("Top 20 Expressed miRNAs") + ylab("Normalized Expression") + xlab("miRNA ID") + 
  theme(panel.grid.major.x = element_line(size = 0.5, linetype = 'solid', colour = "grey"), 
        panel.grid.major.y = element_blank(),
        axis.text.x = element_text(angle = 90, hjust = 1))

ggplot(norm.expr.melt, aes(x=reorder(miR.ID, norm.expr, FUN="median"), y=norm.expr, fill=Race)) + geom_boxplot(pos="dodge", outlier.size=0.5) + 
  ggtitle("Top 20 Expressed miRNAs") + ylab("Normalized Expression") + xlab("miRNA ID") + 
  theme(panel.grid.major.y = element_line(size = 0.5, linetype = 'solid', colour = "grey"), 
        panel.grid.major.x = element_blank()) + coord_flip()

Heatmap of most highly expressed miRNAs

library(pheatmap)
ruvr2 <- ruvr_sets[[2]]
norm.expr.matr <- norm.expr.matr[order(rowMeans(norm.expr.matr), decreasing=TRUE),]
Demographics <- pData(ruvr2)[,c("Age","Race", "Gender")]
annotations <- as.data.frame(Demographics)
rownames(annotations) <- rownames(pData(ruvr2))
pheatmap(norm.expr.matr[0:50,], cluster_rows=FALSE, show_rownames=TRUE, cluster_cols=TRUE, annotation_col=annotations, main="Top 50 Expressed miRNAs")

DE analysis with EdgeR

# based on Ryan's code, pass RUVSeq-corrected data to edgeR
ruvr2 <- ruvr_sets[[2]]
design <- model.matrix(~ 0 + Race + Gender + Age + W_1 + W_2, pData(ruvr2))
dge <- estimateDisp(dge, design = design, tagwise = TRUE, robust = TRUE)
fit <- glmFit(dge, design)

Core miRNAs stable across individuals (Fig2D)

ruvr_norm_matrix <- assayData(ruvr2)$normalizedCounts
#ruvr_norm_matrix <- counts
median <- apply(ruvr_norm_matrix, 1, median)
median_expr_df <- as.data.frame(median)
median_expr_df$qcd <- apply(ruvr_norm_matrix, 1, FUN=function(x){ q <- quantile(x, c(0.25, 0.75)); diff(q)/sum(q) })  # quartile coeff. of dispersion
median_expr_df$cv <- apply(ruvr_norm_matrix, 1, FUN=function(x){ 100*sd(x)/mean(x) })  # coefficient of variation
# log-2 transform the medians
median_expr_df$log2_median <- log(median_expr_df$median, 2)
# plot all miRNAs
ggplot(median_expr_df, aes(x=qcd, y=log2_median)) + geom_point() + ylab('Log2 Median Expression') + xlab('Quartile Coefficient of Dispersion') + ggtitle("miRNA QCD in Plasma")

ggplot(median_expr_df, aes(x=cv, y=log2_median)) + geom_point() + ylab('Log2 Median Expression') + xlab('Coefficient of Variation') + ggtitle("miRNA %CV in Plasma")

# plot only miRNAs expressed in 90% of samples
median_expr_df$sample_fraction = (rowSums(ruvr_norm_matrix > 0) / length(rownames(samples)))
sample_fraction_threshold = 0.9
median_expr_df_filt = median_expr_df[median_expr_df$sample_fraction >= sample_fraction_threshold, ]
ggplot(median_expr_df_filt, aes(x=qcd, y=log2_median)) + geom_point() + ylab('Log2 Median Expression') + xlab('Quartile Coefficient of Dispersion') + ggtitle("miRNA QCD in Plasma")

ggplot(median_expr_df_filt, aes(x=cv, y=log2_median)) + geom_point() + ylab('Log2 Median Expression') + xlab('Coefficient of Variation') + ggtitle("miRNA %CV in Plasma")

Contrast gender (Fig3B)

lrt_gender <- glmLRT(fit, coef = "GenderMale")
alpha_level <- 0.05
results_gender <- data.frame(topTags(lrt_gender, n=Inf, sort.by="PValue", p.value=alpha_level))
sig_miRs_gender = list()
logFC_threshold <- 1
sig_miRs_gender[[1]] <- rownames(results_gender[results_gender$PValue < alpha_level & results_gender$FDR < alpha_level & abs(results_gender$logFC) > logFC_threshold,])  # filter by p-value and logFC
sig_miRs_gender[[2]] <- rownames(results_gender[results_gender$PValue < alpha_level & results_gender$FDR < alpha_level,])
# heatmaps of significant DE miRNAs
Gender <- pData(ruvr2)[,c("Gender")]
annotations_gender <- as.data.frame(Gender)
rownames(annotations_gender) <- rownames(pData(ruvr2))
for(miR_list in sig_miRs_gender) {
  pheatmap(norm.expr.matr[miR_list,], cluster_rows=TRUE, show_rownames=TRUE, cluster_cols=TRUE, annotation_col=annotations_gender, main="Differentially Expressed miRNAs")
}

# TODO: boxplots in one figure
sig_miRs_gender_table <- norm.expr.matr[sig_miRs_gender[[1]],]
sig_miRs_gender_table <- melt(sig_miRs_gender_table)
colnames(sig_miRs_gender_table) <- c("miRNA", "Sample", "Expression")
annotations_gender$Sample <- rownames(annotations_gender)
sig_miRs_gender_table <- merge(sig_miRs_gender_table, annotations_gender, by="Sample")
ggplot(data = sig_miRs_gender_table, aes(x=miRNA, y=Expression)) + geom_boxplot(aes(fill=Gender)) + 
  facet_wrap(~ miRNA, scales="free") + 
  ggtitle("Differentially Expressed miRNAs in Plasma by Gender")

Contrast race (Fig3B)

lrt_race <- glmLRT(fit, contrast=makeContrasts(asn=(RaceAsian-(RaceAfrican_American+RaceWhite)/2),
                                          afr=(RaceAfrican_American-(RaceAsian+RaceWhite)/2),
                                          wht=(RaceWhite-(RaceAfrican_American+RaceAsian)/2), levels=design))
results_race <- data.frame(topTags(lrt_race, n=Inf, sort.by="PValue", p.value=alpha_level))
sig_miRs_race_list = list()
sig_miRs_race_list[[1]] <- rownames(results_race[results_race$PValue < alpha_level & results_race$FDR < alpha_level & abs(results_race$logFC.asn) > logFC_threshold | abs(results_race$logFC.afr) > logFC_threshold | abs(results_race$logFC.wht) > logFC_threshold,])  # filter by p-value and logFC
sig_miRs_race_list[[2]] <- rownames(results_race[results_race$PValue < alpha_level & results_race$FDR < alpha_level,])
# heatmaps of significant DE miRNAs
Race <- pData(ruvr2)[,c("Race")]
annotations_race <- as.data.frame(Race)
rownames(annotations_race) <- rownames(pData(ruvr2))
for(miR_list in sig_miRs_race_list) {
  pheatmap(norm.expr.matr[miR_list,], cluster_rows=TRUE, show_rownames=TRUE, cluster_cols=TRUE, annotation_col=annotations_race, main="Differentially Expressed miRNAs")
}

# boxplots of significant DE miRNAs
sig_miRs_race_table <- norm.expr.matr[sig_miRs_race_list[[1]],]
sig_miRs_race_table <- melt(sig_miRs_race_table)
colnames(sig_miRs_race_table) <- c("miRNA", "Sample", "Expression")
annotations_race$Sample <- rownames(annotations_race)
sig_miRs_race_table <- merge(sig_miRs_race_table, annotations_race, by="Sample")
sig_miRs_race_table <- sig_miRs_race_table[sig_miRs_race_table$Race != "Multiracial" & sig_miRs_race_table$Race != "Pacific_Islander",] # Pacific_Islander and Multiracial have too small sample sizes
ggplot(data = sig_miRs_race_table, aes(x=miRNA, y=Expression)) + geom_boxplot(aes(fill=Race)) + 
  facet_wrap(~ miRNA, scales="free") + 
  ggtitle("Differentially Expressed miRNAs in Plasma by Race")

Create a PCA plot showing Age x Gender

plotPCA(
            reads_sce,
            colour_by = "Age",
            shape_by = "Gender",
            size_by = "total_features",
            exprs_values = "RUVr k= 2"
) + ggtitle("RUVSeq-Normalized Expression (k=2)") 

non-plotting arguments like 'exprs_values' should go in 'run_args'